knitr::opts_chunk$set(warning=FALSE)
knitr::opts_chunk$set(cache = TRUE)For the tutorial below, we are going to use the following R packages:
dplyr for dataframe manipulationggplot2 and ggpubr for plottingGEOquery for handling data stored on NCBI’s Gene Expression Omnibus (GEO)readr for loading large dataframeslibrary(dplyr)
library(ggplot2)
library(ggpubr)
library(GEOquery)
library(readr)Model parameters are stored in Supplementary file 3, but can also be loaded from the following github URL:
https://raw.githubusercontent.com/pipekorsi/MepiClock/main/resources/model_coefficients.csv
df_model_coeff <- read.csv("https://raw.githubusercontent.com/pipekorsi/MepiClock/main/resources/model_coefficients.csv", sep=";", stringsAsFactors = FALSE, check.names = FALSE)The first few raws of the dataframe show its basic structure. Coefficient values equal to zero represent CpG probes not included in the covariates of the given model.
df_model_coeff %>% head()In this tutorial we are going to use one of the methylation datasets included in the test data used for accessing model performance. Its GEO accession number is GSE152026 and in the following we are going to use the GEOquery package to collect relevant information about the study.
gse <- getGEO("GSE152026")We can check the details of the study with the following code:
gse[["GSE152026_series_matrix.txt.gz"]]@experimentData@title## [1] "Blood DNA methylation profiles from first episode psychosis patients and controls I"gse[["GSE152026_series_matrix.txt.gz"]]@experimentData@abstract## [1] "Genome-wide patterns of DNA methylation were quantified using the Illumina Infinium HumanMethylationEPIC BeadChip (“EPIC array”) in DNA samples isolated from blood for schizophrenia cases, first episode psychosis patients and controls. These samples were profiled as part of a wider study where they were meta-analysed with other cohorts."The samples in this dataset are either from patients suffering from schizophrenia or control cases. In our analysis, only the control cases will be needed. We will also need the chronological age of the patients to be able to test model performance. Here we extract age and sample name for control cases and store them in a dataframe:
keep_samples_idx <- grepl("Control", gse[["GSE152026_series_matrix.txt.gz"]][["characteristics_ch1"]], fixed = TRUE)
cc <- strsplit(as.character(gse[["GSE152026_series_matrix.txt.gz"]][["characteristics_ch1.2"]]), ': ', fixed = TRUE)
age <- as.numeric(unlist(cc)[2*(1:length(as.character(gse[["GSE152026_series_matrix.txt.gz"]][["characteristics_ch1.2"]])))])
cc <- strsplit(as.character(gse[["GSE152026_series_matrix.txt.gz"]][["title"]]), ' ', fixed = TRUE)
sample_names <- unlist(cc)[2*(1:length(as.character(gse[["GSE152026_series_matrix.txt.gz"]][["title"]])))-1]
sample_metadata_df <- data.frame(sampleName = sample_names[keep_samples_idx],
age = age[keep_samples_idx], stringsAsFactors = FALSE)
sample_metadata_df %>% head()The number of remaining samples:
dim(sample_metadata_df)[1]## [1] 521Apparently, the Illumina Infinium HumanMethylationEPIC BeadChip was used to generate methylation data, thus we expect to have around 850K measurement points for each sample. Unfortunately, there is no “Series Matrix” available on GEO for this dataset, which usually containes the relevant data in a concise manner, thus the raw data needs to be downloaded and preprocessed.
We downloaded the GSE152026_EUGEI_processed_signals.csv.gz file from the GEO ftp site on a Linux system with:
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE152nnn/GSE152026/suppl/GSE152026_EUGEI_processed_signals.csv.gz
Checking the first few rows and columns of the downloaded datafile reveals that both the \(\beta\)-values and the detection p-values are recorded and samples are referred to by their sentrix position and chip ID:
as.character(system("zcat data/GSE152026/GSE152026_EUGEI_processed_signals.csv.gz | head -2 | cut -d',' -f1,2,3,4,5"))"","200863670076_R01C01","200863670076_R01C01_Detection_Pval","200863670076_R02C01","200863670076_R02C01_Detection_Pval"
"cg00000029",0.37937352230886,1.62762280961683e-26,0.444926666863015,1.47456418339601e-85Checking the number of lines is in agreement with our previous knowledge that the samples were processed on the EPIC array:
system("zcat data/GSE152026/GSE152026_EUGEI_processed_signals.csv.gz | wc -l")860964Given that the epigenetic clock models only contain a couple hundred CpG sites, it is computationally more efficient to only load a dataset that was previously filtered for these probes and sample columns:
probes_to_keep <- df_model_coeff$probeID[df_model_coeff$probeID != "(Intercept)"]
samples_to_keep <- sample_metadata_df$sampleName
f <- function(x, pos) x %>% select(c("X1", all_of(samples_to_keep))) %>% mutate(probeID = X1) %>% filter(probeID %in% probes_to_keep) %>% select(-X1)
beta_df <- read_csv_chunked("data/GSE152026/GSE152026_EUGEI_processed_signals.csv.gz",
DataFrameCallback$new(f),
chunk_size = 10000,
quote = "\"",
progress = TRUE
)